
colors = {'#000000','#0072BD','#A2142F'};
colors={[0 0.4470 0.7410],[0.8500 0.3250 0.0980],[0.9290 0.6940 0.1250],[0.4940 0.1840 0.5560],[0.4660 0.6740 0.1880],[0.3010 0.7450 0.9330],[0.6350 0.0780 0.1840]};
shapes = {'s','^','o'};
% Q2
% 22 + 0,0,-4,1/2, '20230417_100VR_30G_10kHz_REMPI_scan_5'
% cyclescut = [59 62]
delay = [ 0    40    80   120   160   200   300];
y2 = [57.2222   36.4000   22.4000   22.7000   17.0000   12.5556    8.3000]-2;
ey2 = [2.5264    1.9183    1.5166    1.5330    1.3115    1.2019    0.9539];


% Q2
% % 22 + 0,0,-4,3/2, '20230629_100VR_30G_10kHz_SR_scan_3'
% delay=[ 0    40    80   120   160   200   300];
% y2=[ 47.4000   34.2000   20.2000   12.0000    7.8000    8.0000    6.6000];
% ey2 = [  3.1048    2.6458    2.1071    1.6248    1.3115    1.2961    1.2490];


% Q1
% 22 + 0,0,-4,1/2, '20230417_100VR_30G_10kHz_REMPI_scan_6'
% cyclescut = []
y1= [25.1667   15.7273   10.6364    9.5833    6.1667    5.0000    3.0909]-2;
ey1 = [1.4860    1.2365    1.0245    0.9317    0.7638    0.7873    0.5750];


% R0
% 22 + 0,0,-4,1/2, '20230417_100VR_30G_10kHz_REMPI_scan_7'
% cyclescut = []
y0 = [27.7000   18.6000   10.6000   11.0000    6.9000    4.6000    3.3000];
ey0 = [1.6882    1.3856    1.0677    1.0770    0.9000    0.7483    0.6403];


% background for R0
% '20230418_100VR_30G_10kHz_REMPI_scan_1'
% cyclescut = [35 71]
bk = [7.1000    4.3000    3.4000    2.5000    1.9000    2.0000    2.0000];
ebk = [0.9000    0.7141    0.6000    0.5385    0.5568    0.5831    0.4690];
y0 = y0-bk;
ey0 = sqrt(ey0.^2+ebk.^2);

Data = {y0, y1, y2};
Error = {ey0, ey1, ey2};
% 
xft = linspace(0,max(delay)*1.01,100);
c0 = [50 0.00001];
figure(90)
hold on
t_half = zeros(1,length(Data));
e_t_half = zeros(1,length(Data));
alpha = zeros(1,length(Data));
e_alpha = zeros(1,length(Data));
for i=1:length(Data)
    t = delay;
    data = Data{i};
    err = Error{i};
    [c,cerr,fit_val] = fit_wrapper(@(c,t) TwoBody(c,t),c0,t',data',err',xft);
    fit_val = TwoBody(c,xft);
    errorbar(t,data,err,shapes{i},'markersize',8,'linewidth',2,'color',colors{i},'displayname',['N',num2str(i)])
    plot(xft, fit_val,'-','Linewidth',1.5,'color',colors{i})
    ct= 1/(c(1)*c(2))*sqrt((cerr(1)/c(1)).^2+(cerr(2)/c(2)).^2);
%     t_half(i) = 1/(c(1)*c(2));
%     e_t_half(i) = ct;
    t_half(i) = c(2);
    e_t_half(i) = cerr(2);
    alpha(i) = c(1);
    e_alpha(i) = cerr(1);
end
t_half 
e_t_half
alpha
e_alpha
legend
box on
set(gca,'fontsize',16)
set(gcf,'color','white')
xlabel('Detection delay (ms)')
ylabel('KRb ion counts per cycle')


function [c,cerr,fit_val] = fit_wrapper(func_hand,c0,x,y,ye,xft)
    if nargin < 5
        ye = 0;
    end
    if nargin < 6
        xft = x;
    end
%     opts = optimoptions('OptimalityTolerance',1e-8);
%     opts = optimset('Display','off');
    opts = optimoptions('lsqcurvefit','OptimalityTolerance',1e-10,'FunctionTolerance',1e-10);
    if ye == 0
        [c,R,J] = nlinfit(x,y,func_hand,c0,opts); % slightly faster
    else
        [c,~,R,~,~,~,J] = lsqcurvefit(func_hand,c0,x,y,[0.0001,0.00001],[5000,3000],opts);
    end
    sum(R.^2);
    ci = nlparci(c, R,'jacobian',J,'alpha',.3173); % 1 sigma
    cerr = abs(ci(:,2)-ci(:,1))/2;
    fit_val = func_hand(c,xft);
end



function y= TwoBody(b,x)
    y=b(1)*(1-x./(b(2)+x));
end


function color = hex2RGB(hex)
    color = sscanf(hex(2:end),'%2x%2x%2x',[1 3])/255;
end